function sT(f,fi)
 global T Ti T0i T0 wi w wi0 C2 C3 vex

   fxy=convect(T0i)+convect(Ti);
%    fxy=convect(Ti);

   fxy = fxy - dif3(T);
%    fxy = C2*fxy -0.05*dif3(T);

   work=T*f+Ti*fi-fxy;
  %%
   if(f<0.6)%the 1st time
       T=Ti;
       Ti=work;
   else %the 2nd time
       Ti=work;
   end
end